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^5 The 3D Airy beam (AiB) is thoroughly explored from a wave-theory point 

g of view. We utilize the exact spectral integral for the AiB to derive local 

I— I ray-based solutions that do not suffer from the limitations of the conventional 

1— I parabolic equation (PE) solution, and are valid far beyond the paraxial 

zone and for longer ranges. The ray topology near the main lobe of the 
^ AiB delineates a hyperbolic umilic diffraction catastrophe, consisting of a 

CO cusped double-layered caustic, but this caustic is deformed in the far range 

^ where the field loses its beam shape. The field in the vicinity of this caustic 

^ is described uniformly by a hyperbolic umilic canonical integral which is 

IL" structured explicitly on the local geometry of the caustic as obtained from 

*^ the initial field distribution. In order to accommodate the finite-energy AiB 

we also modify the canonical integral by adding a complex loss parameter. 
The canonical integral is calculated using a series expansion and the results 
are used to identify the validity zone of the conventional PE solution. The 
analysis is performed within the framework of the non-dispersive AiB where 
the aperture field is scaled with frequency such that the ray skeleton is 
frequency-independent. This scaling enables an extension of the theory to the 
ultra wide band (UWB) regime and ensures that the pulsed field propagates 
along the curved beam trajectory without dispersion, as will be demonstrated 
in a subsequent publication. 
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1. Introduction 

Airy beams (AiB) are exact solutions of the time-harmonic parabolic wave equation (PE) in 
free space, that are generated by an Airy function initial condition in an extended aperture. 
The AiB's have attracted a considerable attention mainly because of their intriguing features, 
the most distinctive one is the propagation along curved trajectories (l}]4]- Other intriguing 
properties are the weak diffraction |2] and self-healing (i.e., regeneration of the solutions if the 
beam is partially obstructed) [s] . Many interesting applications that utilize these properties 
have been suggested, e.g. l6|-[To]. 



In a previous paper 11 , we provided a cogent wave-based description to the 2D- AiB's 



that clearly explains all these features. Using uniform asymptotic evaluation of Kirchhoff 's 



radiation integral 12 and uniform geometrical optics (UGO) [13 ,14 , it has been shown that 
the AiB is not generated by the main lobe of the Airy function in the aperture, but it is 
in fact a caustic of rays that emerge sideways from the side lobes of the extended aperture 
distribution. This implies that the AiB cannot be regarded as a "beam field" in the sense 
that the evolution of the main lobe along the curved propagation trajectory is not described 
by a local wave dynamics. It has also been established that ray theory in conjunction with 
the UGO provides an accurate solution for the AiB, which fully agrees with the conventional 
PE solution of []1,'2] when the latter is valid. The UGO solution is accurate even at large 
observation ranges where the PE solution fails due to the curving beam axis. It has been 
noted that these ray-based tools can be used to construct other curving beam solutions in 
uniform and non-uniform media by back-projecting the rays from the caustic to the aperture 



and synthesizing the aperture distribution that focuses onto that caustic 11 . This concept 
has been applied recently in [Tsl. 



Using the ray description of |11] we have also constructed in 16 a class of non- dispersive 
Airy pulsed beams (AiPB). These solutions are obtained by synthesizing the aperture field 
distribution such that the ray skeleton is the same for all frequencies. This scaling of the 



aperture distribution, unlike previous suggestions 17 , ensures that all frequency components 
focus on the same curved caustic, thus giving rise to a pulse that propagates along the curved 
trajectory with no dispersion. A closed-form time-domain expression for the AiPB has also 
been derived in [l6] using the spectral theory of transients (STT) [18 -20 . 

The main goal of this paper is to clarify the wave mechanism of the 3D-AiB by deriving 
uniform ray-based solutions that reveal the local features of the field. Furthermore, these 
solutions rely on the exact ray and caustic topology of the AiB and therefore do not suffer 
from the limitations of the conventional PE solution. We note that due to the curving of 
the beam trajectory, the PE solution is valid only for short propagation distances and in 
particular it does not predict the termination of the AiB at a certain range, as will be 
discussed here. In order to meet this goal we first analyze the ray topology of the 3D- 
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AiB (Sec. [3]) and then use this analysis to derive uniform solutions that are structured 
explicitly upon the local geometry of the caustic (Sees. 4-5). As will be explained below, 
we emphasize here local uniform solutions rather than globally uniform solutions since the 
former yields explicit solutions that explain all the parameters and features of the field. 
The analysis is performed within the framework of the 3D non- dispersive AiB, defined in 
Sec. [2] such that the ray skeleton is frequency-independent. The non-dispersive formulation 
not only simplifies the ray analysis by making it more transparent, but it also enables the 
construction of the non-dispersive AiPB as noted above. This subject will be explored in a 
subsequent publication f2T], where the results of the present work will be used to synthesize 
3D-AiPB and to derive closed-form time-domain solutions. 

In Sec. |3]we study the exact ray and caustic topology of the 3D-AiB, concentrating on the 
region near the curved beam axis and deferring the discussion on the global caustic topology 
to Sec. [6j The analytical details are summarized separately in Appendix |Aj The topology 
near the beam axis has many intricate features and it is considerably more complicated than 
the approximate ray topology of the 3D PE solution in [7j or the 2D case in 11 . Specifically, 
we show that the main lobe of the AiB is described by a cusped double-layered caustic which 



is of the hyperbolic umbilic type 22 26] (see Figs. |4]-[5]). We also identify the range where 
the isolated hyperbolic umbilic caustic that describes the AiB field near its propagation axis, 
evolves into the more complicated parabolic umbilic caustic considered in Sec. [6} so that the 
AiB loses its form. It should be noted that these caustic evolutions (unfoldings) are expected 
from catastrophe theory principles (see e.g.. Fig. 4.4 and Fig. 7.3.(3-5) in |24], respectively), 
yet, in this paper we provide exphcit expressions for this overall topology and also use it to 
construct the field solution, as described below. 

The field calculations are considered next in Sees. |4]-[5| We start from the exact spec- 
tral diffraction integral for the AiB, and then discuss its asymptotic evaluation in terms of 
geometrical optics (GO). The field near the beam axis, where the ray topology is of the hy- 
perbolic umbilic type, is considered separately in Sec. |5] within the framework of diffraction 
catastrophe theory [24|[27|[28] , discussed below. 

Catastrophe theory [24||27||28] classifies the caustic topologies in stable polynomial forms, 
referred to as elementary catastrophes. This also provides a systematic framework for cal- 
culating the field near caustics with arbitrary topologies. The general methodology is to 
map the original integral representation of the field, expressed either spatially or spectrally 
(e.g., Kirchhoff or plane- wave integrals, respectively), onto the diffraction catastrophe inte- 
gral corresponding to the elementary catastrophe that describes the caustic. This integral, 
also referred to as canonical integral, constitutes a standard special function 29, Sec. 36.2]. 

The mapping of the diffraction integral into the canonical catastrophe integral may be 
applied globally via a systematic method first introduced in 30, Sec. 5]. This method ex- 
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tends to multi-dimensions the method of Chester et al. fr2] (see also 31,32]) for uniform 



asymptotic evaluation of ID integrals with several nearby stationary points and singulari- 



ties 31 . The global mapping approach, while very elegant and systematic, leads to implicit 
representations where the field is described in a generalized coordinate system that has to 
be calculated numerically. Alternatively, this mapping may be performed locally around any 
given point in configurational space, leading to explicit field representations where the gen- 
eralized coordinates are the local Cartesian coordinates. In this paper we utilize the local 
mapping approach since our goal is to explore the wave mechanism of the 3D-AiB by de- 
riving explicit field solutions where all the parameters are directly related to the aperture 



field distribution. We start in Sec. 5. A by mapping the exact AiB spectral integral in the 



vicinity of the curved beam axis, ending up with a hyperbolic umbilic diffraction integral 



for the field. Next in Sec. 5.B we extend this integral to the case of finite- energy AiB's by 



adding a complex loss parameter into the phase of this integral. In Sec. |5.C we calculate this 



integral by extending the closed-form series expansion of 29, Sec. 36.8] to finite energy AiB's 



Finally, in Sec. 5.D we use this solution in order to check the validity of the conventional PE 
solution, and demonstrate that the latter suffers from a translational error due to an error 
in the beam axis under the PE approximation. 

2. Non Dispersive Airy Beams 

In a previous paper |11| we have dealt with the ray interpretation to the 2D-AiB and demon- 
strated that its curved propagation trajectory is a caustic of rays that emerge sideways from 
an extended aperture. Based on this concept we have extended the 2D-AiB to the ultra wide 
band (UWB) frequency domain and to the short-pulse time domain, by employing an initial 
field distribution that gives rise to a frequency-independent ray skeleton, thus ensuring that 
all frequency components propagate along the same trajectories and focus onto the same 



caustic 16 



Extending this concept to the 3D case, the non-dispersive AiB field in the half-space z > 
is generated by the aperture field distribution in the z = plane 

f/o(x,y;w) =exp[fcKx + «,y)]Ai(/3;i/=^A;2/3a;)Ai(/3-i/3^2/3^), (1) 

where (x, y) are the transversal coordinates, Ai is the Airy function, and a harmonic time- 
dependence exp(— iwt) is assumed and suppressed throughout, k = u/c is the wave-number 
and the wave-speed c is assumed uniform. The exponential window in Eq. ([T]) is added 
in order to render the energy of this distribution finite. The parameters Px,f3y, cxxiC^y are 
frequency independent scaling factors, with Px,y controlling the widths of the field distribution 
along the x and y axes, respectively, and ax,y controlling the decay rates along these axes 
(see Fig.g. 



4 



It should be noted that the conventional expression for the AiB initial condition of Eq. ([T]) 
is Uo{x,y) = Ai{x/x„)Ai{y/yo) exp{axX + ayy) fs", Eq. (2)]. The scaling parameters in this 
expression are related to those in Eq. ([T]) via Xq = pl^^k"^/^, y^ = P^^k''^/^ and a^.y = kax,y 
We prefer the frequency independent (5x,y parameterization since it renders the ray skeleton 
and the propagation trajectory frequency independent (see Eq. ( [Io| and Sec. [s]). 

The radiating field at r = (x, z) due to the distribution in Eq. ([T| has a closed- form 
solution within the parabolic wave equation (PE) which can be expressed as [s] 



U (x, y, z) = F{x, z; u, ax)F{y, z; u, /3y, ay) exp{ikz), (2) 

where F is given by 

F{p,z-u,f3,a) = Ai[{kf3y/''{p/f3 - {z/2f3f + iaz/f3)] 
exp [ik{pz/2(3 - z^/l2(3'^ + a^z/2)\ 
exp [A;a(p-zV2/3)]- (3) 

The beam trajectory is defined by the vanishing of the arguments of the Airy functions in 
Eqs. (|2|-(|3|), and is given bjQ 

xll3x - {zl2(3xf = 0, (4) 
y/Py - {z/2Pyf = 0. (5) 

Note that the Px,y parameterization implies that the beam trajectory is frequency indepen- 
dent. 

For later use we introduce the coordinate system {x,y,z) such that the beam trajectory 
is contained in the {x,z) plane (see Fig. [T]). This system is obtained by rotating the {x,y) 
system at the angle ip = ta.n~^ {f3x / f3y) about the z-axis, viz 

{x,yy = R^{x,yY, (6) 

where the superscript T defines the transpose of a matrix and is the rotation matrix 

R^^l""" """V (7) 

\— sm(y9 cos if J 

For simplicity, but without loss of generality, we shall consider the symmetrical case 

/3. = /3, = /3 = ^^2, (8) 
Oix = oiy = a = a/ "\/2, (9) 



^ Actually, the peak trajectory is described by replacing the on the right hand side of Eqs. (|4])-([5]) by the 
asymptotically small parameter (fc/3)~^/'^a']^ where a'^ ~ —1.0188 is the first zero of Ai'(a;) 29 Table 9.9.1]. 
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such that V9 = 45° (see Fig. [T]). The parameters in the (x, y) system are defined conveniently 
as P and a, and the beam trajectory (the "beam axis") is given by 

x/(3 - {z/2(3Y = 0. (10) 

The general case ^ f3y can be reduced to the symmetrical case by the scaling x — )■ 
^[{/^x + l^y)^^'^ / f^xY^^ , y y[{Px + l^yY^^ / f^yV^^ ^"^^ applying the exponential windows in the 
new coordinate system. 

The solution in Eqs. (|2])-([3]) is derived within the PE (paraxial) approximation and is 
therefore inaccurate at large distances where the AiB propagation trajectory curves. The 
ray solution of Sees. |4]-[5} on the other hand, remains accurate even at large distances since 
it is structured upon the ray skeleton of the exact wave equation. The ray solution also 
provides a clear physical interpretation for the AiB and a tool for constructing other curving 
beam solutions in uniform and non-uniform media. 

3. Ray System and Caustic Topology 

3. A. Ray System 

The properties of the AiB are fully described by the ray system that emerges from the 
aperture plane. We therefore start with the ray-system topology. Position r = (x,y,zY 
along a ray can be expressed as 

r = r +S0-, (11) 

where r' = (x',?/',0)^ is the ray initiation point in the z = plane, a is the distance along 
the ray, and s defines the ray direction viz 

§ - (^7 ^, C)^ - (cos6'x., cos 6y, cosO^y, (12) 
with being the ray angles with respect to the x,y,z axes, respectively, and 

c = {i-e-vY'- (13) 

In order to find s for a given initiation point (x', y') we substitute in Eq. ([!]) the asymptotic 
form of the Airy function, 

Ai(a;) ~ (-XTr^)-!/^ sin [^{-x)^^^ + x < -1, (14) 

and recast the initial field in the form 

4 

Uo{x',y') ^ ^Aosexp[ikipos\, (15) 

s=l 
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where 



^^1/6 ^ka{x'+y') 



,-iTTMos/2 



AiTk^^ i-x'y/'^i-y'y/'^ 

with the ± signs in Eq. (16) taken as ++, H — , — h, , for s 

for s 



(16) 
(17) 



Eq. (15), Mo, 



1, Mos = 1 for s = 2, 3 and Mos 



1,2,3,4, respectively. In 
2 for s = 4. We assumed here 



that a is small so that the decay term exp[ka{x' + y')] is slowly varying on a wavelength 
scale and can be included in the amplitude A^g rather than in the phase term ipQg. This 



assumption is justified in Sec. 4.B 



One finds that each point {x',y') in the aperture plane, with x',y' < 0, emits four rays, 
corresponding to s = 1, .., 4 in Eq. (15). The ray directions are obtained by local matching 
to the phase gradient viz 

[e., Vs] = [d.', dy,]As{x',y') = [± (-x7/3)l/^ ±{-y'/(3Y/'], (18) 



where the "±" signs are taken as explained in Eq. (16). Note that the ray directions are 
frequency independent if /3 is frequency independent. The directions of the four ray species 
projected onto the z = plane are depicted in Fig. [2j Finally, we note that the result in 
Eq. (18) has been derived asymptotically via Eq. ( [I4) ), but as shown in Appendix A.l, it is 
also valid exactly. 



The ray trajectories in Eq. (11 ) are parameterized in terms of the ray coordinates (^, t], a) 
where the ray direction {^,1]) tags the initiation point {x',y') via Eq. ( [l8| . The topology of 
the ray system is described by the Jacobian of the mapping from the (x, y, z) coordinates to 
the ray coordinates (^,77,0"), defined by 

d{x,y,z) 



J 



(19) 



where = [d^x, d^y, d^zY and || || denotes the determinant. An explicit expression for J 
is obtained by noting from Eq. (18) that the exit point coordinates are defined by the ray 



coordinates via {x\y') 
in Eq. (19) we obtain 



J 



Substituting in Eq. (11) and calculating the derivatives 

7/ a-2/3r/ . (20) 

0. This condition defines the envelope of the ray system, 
also termed caustic surface or catastrophe 24 ,28 . The structure of the caustic near the beam 
axis is discussed next in Sec. |3.B , with further details in Appendix |X} The complete caustic 
topology including also regions away from the beam axis is presented in Sec. |6] 



The mapping is singular when J 
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3.B. Caustic Topology in the Vicinity of the Beam Axis 
3.B.I. Caustic topology within the PE model 



The PE approximation for the beam axis in Eq. ( 10 ) is vahd only for moderate ranges before 



the axis curves too much. In this range one may approximate in Eq. (13) C ^ 1 — ^(.^^ + 77 



Substituting into Eq. (20) one obtains a closed form solution for the caustic equation ^7 = 



see Appendix A.3). Figure 3 depicts the caustic structure within this approximation. It 



consists of two curved sheets joined at an edge that delineates the beam axis of Eq. (10). 
Fig. |4] shows a cross sectional cut of the caustic in two constant-z planes, where it is seen to 
form a right-angled corner. Figure |4] also depicts the exact caustics to be discussed next. 

3.B.2. Exact caustic topology along the beam axis 



Employing exact ray theory (see Appendix A. 2) reveals further details of the caustic, which 



are also evident from the canonical catastrophe integral to be discussed in Sec. |5} Here we 
discuss the structure of the caustic near the beam axis; the complete global structure beyond 
this zone will be discussed separately in Sec. [6| 

The exact caustic is shown in Fig. |5} It consists of two surfaces, referred to as "caustic 1" 
and "caustic 2," each consisting of two sheets joined at an edge. Caustics 1 and 2 are very 
close and hardly distinguishable on the figure scale. The caustic within the PE approximation 
is a degenerate case in which caustics 1 and 2 coalesce. Figure |4] zooms on cross-sectional cuts 
of the caustics at two constant-2; planes, showing that caustic 1 is a cuspoid, whereas caustic 2 
is smooth, and that the distance between these caustics increases with z. This topology is 
identified as a hyperbolic umbilic diffraction catastrophe (see e.g.. Fig. 4.4 in [24]) whose focal 
plane is located at the z = plane where the two surfaces coalesce along the semi-infinite 
lines {x = 0,y < 0} and {y = 0,x < 0} and form a right-angled corner at the origin. 

The explosion of the PE caustic into a hyperbolic umbilic caustic at z > can also be 



expected from catastrophe stability theorems 27,28 as follows from the fact that the right- 



angled corner at the cross-section of the PE caustic is not a stable form (the only stable 
forms in 2D are the fold and the cusp; see e.g.. Table 3.1 in [24]). Therefore, the corner will 
not survive the perturbation due to exact propagatioi][^ Furthermore, the simplest stable 
singularity in 3D which contains a corner at a particular cross-section (termed focal plane) 
is the hyperbolic umbilic catastrophe (see also |33] where a similar argument has been used). 

Figure |4] also shows that the deviation of the PE caustic from the exact one increases 
with z, thus limiting the range of validity for the PE solution in Eqs. ([2])-([3]). In Fig. |6|we 
concentrate on the difference between the edges of the exact caustics 1 and 2 (Eqs. ( A.16[ )- 



(A.19)) and of the PE caustic (Eq. (10)) since this edge delineates the propagation trajectory 



The PE model is a special case where the effect of the propagation is only a translation (see Eqs. (A. 23 ) 



(A.24|) and therefore the caustic retains its unstable shape. 
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of the AiB. 

Figures [7] and [8] describe the ray formation of caustics 1 and 2, respectively, by depicting 
the rays emitted from points in the aperture plane and their point of tangency at the caustics. 
Actually, the figures show traces of the exit points in the aperture (thick lines), and the traces 
on the caustics of the corresponding points of tangency (thin lines with the arrows indicating 
the direction of the traces). To emphasize the symmetry, traces with y' ^ appear as full 
or dashed lines, respectively. Lines of constant y' and of constant x' are plotted in magenta 
and blue, respectively. 

As discussed in Eq. (18) and Fig. [2| each point in the aperture emits four rays, belonging 



to species s = 1, ...4. As seen in these figures, caustic 1 is formed by ray species 1,2,3, whereas 
caustic 2 is formed by species 1 only. Specifically, referring to Fig. [7} the rays of species 2 
and 3 touch only caustic 1 (red points), whereas the rays of species 1 penetrate through 
caustic 1 (through the white circle), then touch caustic 2 (green point; see also Fig. [s]), 
penetrate again through caustic 1 (white circle) and finally touch caustic 1 (red point). In 
the limit — )■ 0, the rays of species 1 tend to the symmetry plane y' = thus giving rise 



to the somewhat counterintuitive situation in Fig. 10 where the ray in the symmetry plane 
penetrates through the edge of caustic 1 and touches caustic 2, but then as it penetrates 
back through the cusped-edge of caustic 1 it is also tangent to the sheets that generate the 
cusp (see the blue solid-line trajectory in Fig. [?]). This situation will be discussed again in 



connection with Figs. M and 10 



4. Exact Spectral Representation of the Field and Geometrical Optics 

4- A. Spectral Representation 



Using the spectral representation of the Airy function 29, Eq. (9.5.1)], the initial field in 
Eq. ([T]) can be expressed as 

^2/3 roo nco 

f/o(x, y-u) = j—yj I d^dvA e--o{€,.)e-«.-+..)/c^ (21) 



OO — OO 



where exp[ia;(^x + riy)/c] is the Fourier kernel, (C,^) are the spectral variables associated 
with (x, y) and 

T-o(e, V) =m + ^«)' + iv + taf]/3c, (22) 
A = {(3/cf/\ (23) 



The spectral integral in Eq. (21) has been expressed in a normalized non-dispersive form |18 
where u appears as a linear parameter in the phase and to is frequency independent. As a 
result (^, Tj) have a pure frequency-independent geometrical interpretation in terms of the 
plane wave angle (see discussion below), and ro is the departure time of the rj) plane wave 
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from the z = plane. We use this normahzation not only because its physical transparency 
but also since it may be extended to the time-domain 



21 



The radiated field, obtained by adding the spectral propagator, is 

^2/3 roo nco 



— oo J — oo 



(27r)2 

= [To{^,r]) +^x + r]y + Cz]/c 



(24) 
(25) 



where ( = yl — — "/^^ is the spectral wave-number in the 2;-direction, chosen with Im C > 
for u > 0. Consequently the integration extends along the real ^, r] axes from — oo to oo such 
that the ^ contour passes above and below the branch points ^ = =Fa/1 — corresponding 
to C 



Equation (24) expresses the field as a spectrum of plane- waves propagating in the direction 



s(^, T]) defined in terms of (^, rj) via Eq. ( 12 ). Note that in Sec. k3l rf) were parameters that 



define the ray in Eqs. ( 11 )-( 12 ), whereas in this section, they are spectral variables that define 



the spectral-plane-wave direction. For simplicity we use the same notations and distinguish 
between them from their context. Actually, the two are related since the ray coordinates of 



Sec. [3] are stationary points in the spectral (^,?7) domain, as explained after Eq. (26) and 
also in Appendix |A.1[ 



As noted after Eq. (21), the spectral integral Eq. (24) has a non- dispersive form where u 



appears as a linear parameter in the phase and r is frequency independent, describing the 
arrival time of the plane wave that passes through r. This special form will be used in [21] 
to evaluate the field of the Airy pulsed beam (AiPB). 



The direct numerical evaluation of the spectral integral in Eq. (24) is difficult since the 



integrand is very oscillatory at high frequencies. It may be evaluated asymptotically using 



the steepest descent method 34 which yields the geometrical optic (GO). Equation (24) will 
be also used in Sec. [5] in order to calculate the AiB within the framework of Catastrophe 
theory. 

4.B. Geometrical Optics 



The main contributions in the integral of Eq. (24) come from the vicinity of the saddle points, 
defined by 



d,T = 0, 



0, 



(26) 



with r given in Eq. (25). We denote the solutions of Eq. (26) as {^r,Vr) where r is the ray- 
index. For observation points far away from the caustic, the saddle points are distinct and 



the field can be expressed in terms of the isolated saddle point contributions 34 

A V- 



U 



-l/3_ 



27T 



g-i(V4K|^ 1-1/2 g 



lUJTr 



(27) 
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with Tr = T a.t {Cr,Vr)- is the determinant of the Hessian matrix 



H 



9|r d^dr^T 



(28) 



evaluated at (C,r,r]r), and fir is determined by the eigenvalues h^r and /lar of H, evaluated at 
i^r,Vr), via 



fir = sgn /lij. + sgn /la 



(29) 



Expressions for {^r,Vr) and H are given Eqs. (A.1)-(A.4) and (A.7)-(A.9), respectively, 



Expressions for Hr, hi 2r and fir are given in Eqs. (A.25)-(A.27) 



As discussed in Eq. (15), a is relatively small and hence we consider first the case a = 0. 



In this case, the saddle points {^r,Vr) corresponding to an observation points r on the lit 



side of the caustic are real and equal to the ray parameters in Eq. (12) of the rays reaching r 



see proof in Appendix A.l ). The ray-index r should not be confused with the species-index 



s introduced in Eq. (15): depending on r, the rays described by the rth saddle points can be 
of any species. 



The result in Eqs. (27)-(29) can therefore be expressed in the GO form (see details in 



Appendix A. 4) 



Ty> — Tq^ ~1~ dr /c, 

j\ g-*('r/4)Atr 



UJ 



-1/3. 



-iMrTT/2 A 
e /In 



j;(o) 



(30) 
(31) 



where is the ray length in Eq. (11) and Jr is the Jacobian in Eq. (19), both evaluated 



at {C,r,Vr)- All squared roots are taken as positive. Mr is the Maslov index that counts the 



number of times ray r has touched the caustic. Tor and Aor in Eqs. (30)-(31) are, respectively, 
the initial delay and amplitude of the rth ray, given by Eq. (16)-(17), with r^r = ipos/c, both 



are taken at {x'j.,yl), the initiation point of the rth ray. Note that the GO amplitude in 
Eq. ( pT| is invalid at the caustics where i7r(c"r) = 0. 

Fig. |9] depicts an illustrative example of the rays reaching a given observation point r. For 
simplicity we consider a point in the vicinity of the beam trajectory where all rays are of 



species 1 (see Fig. 2(a) ). A point on the lit side of the caustic is reached by four rays, tagged 
by the index r = 1,2,3,4. The Maslov index of Eq. (31) in this example is as follows: Ray 
r = 4 has not touched the caustic hence M4 = 0; Rays r = 2, 3 have penetrated through 
"caustic 1," then touched "caustic 2", and then penetrated again through "caustic 1" (see 
also ray 1 in Fig. [T]), hence = 1; Ray r = 1 is has penetrated through "caustic 1", 
touched "caustic 2" and then touched the cusp at the edge of "caustic 1" while penetrating 
it (cf. Fig. 10), so Ml = 2. In addition to the four rays described above, the equation for the 
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rays has three more solutions representing rays that emanate from the aperture at distant 
points where the field is very weak and the ray angles are very shallow with respect to the 
z = plane. These rays will not be considered therefore in the analysis. 



Finally, we note that the GO field in Eq. (27) has been derived via an asymptotic analysis 
of the exact spectral representation, but for a = it reduces exactly to the asymptotic 



initial conditions in Eq. (15) that were derived from the asymptotic approximation of the Ai 



function in Eq. (14). 



If a 7^ then the saddle points become complex. The saddle point analysis above still 
applies, but the physical interpretation is less transparent. However, since typically a is very 



small, one may expend and Hr of Eq. (27) into Taylor series about the real saddle points 
{^r, Vr) for <y = 0. Keeping terms to leading order in a one obtains the result in Eq. (27) with 



Eqs. (30)-(31) except that here Agr — )■ AorC (cf. Eq. (17); see details in Appendix A. 4) 



4-0. The field Near the Caustic 



Near the caustic, the ray-system Jacobian vanishes and the GO amplitude in Eq. (31) ex- 



plodes. From a spectral perspective, the relevant saddle points come close together and the 



second order asymptotic analysis in Eqs. (27)-(29) is no longer valid. One should resort 



instead, to higher order asymptotics that may accommodate two or more nearby saddle 
points. The simplest situation occurs near the smooth parts of the caustics which are fold 
catastrophe where two saddle points coalesce. In these regions, the field may be described 



uniformly via the uniform geometrical optics (UGO) [13, 14 where the field is expressed 
explicitly in terms of the ray-parameters (i.e., the phase and the Jacobian), but these pa- 
rameters are combined together so that the solution remains regular near and on the caustic 
while reducing to the conventional GO expressions further away from it. The simple UGO 
is not valid near the beam trajectory where the caustic topology is a hyperbolic umbilic as 
discussed in Sec. |3.B.2| and Figs. |4]-[8j The field in that region is considered next in Sec. [5] 



5. The Field Near the Propagation Trajectory 

Near the beam trajectory, the caustic topology is that of a hyperbolic umbilic, hence, as ex- 
plained in the introduction, we refer to the general framework of catastrophe theory in order 
to calculate the field. In this approach, the field is represented by the canonical catastrophe 
integral that can be written as 

Co 



U{x,y,z) 



drj e 



(32) 



where {x, y, z) are generalized coordinates in the vicinity of the caustic such that z is a 
coordinate along the caustic and (x, y) are transversal coordinates, and (^, ff) are the spectral 
variables associated with {x,y). $ has a canonical polynomial form for each elementary 



12 



catastrophe (see e.g., Eq. (42) for the hyperbohc umbihc). Near the caustic, the integral in 



Eq. (32) can be evaluated by standard numerical techniques such as series expansion 29 



Sec. 36.8] (see also discussion in Sec. 5.C), whereas further away from the caustic, it reduces 
asymptotically via saddle point analysis to the GO field. The spatial and spectral coordinates 



in Eq. (32) are typically normalized to be dimensionless and Cq is a normalization constant 
such that t7 = 1 at [x, y, z) 



0. 



In the theory of diffraction catastrophes 23 , 30 , the exact integral representation of the 



field, expressed either spatially (i.e., via a Kirchhoff-type integral) or spectrally, is mapped 



asymptotically via a method outlined in 30, Sec. 5] into the canonical integral in Eq. (32), 
thus leading to a globally uniform solution of the fieldj^ 

As discussed in the introduction, in the globally uniform approach, the generalized coor- 
dinates in $ are related implicitly to the geometrical coordinates near the caustic. Our goal, 
however, is to obtain an explicit mapping such that these generalized coordinates are directly 
related to the local Cartesian coordinates along the beam trajectory as shown in Fig. [10} 
and also to the parameters of the aperture field distribution. Such explicit expressions are 



obtained by the locally uniform approach of Sec. 5. A where we concentrate on the region 



near the caustic. In this approach, we transform the exact spectral integral of Eq. (24) to the 



local Cartesian coordinate system in Fig. 10 , and then approximate the spectrum paraxially 



around the caustic direction. This leads directly to the spectral integral in Eq. (32) where the 



phase is recognized as hyperbolic umbilic. This integral describes the field in the transition 
layer around the caustic, but its validity range is much larger and in particular it reduces to 
the GO field away from the caustic. As discussed in the introduction, this direct approach 
also provides a convenient tool for the synthesis of other types of beams. 

As was done in Sec. |3| we consider first in Sec. 5. A the case where a = where the rays 



are real. The modifications for the finite energy AiB where a 7^ are done in Sec. 5.B In 



Sees. 5.C 5.D we discuss the numerical evaluation of the integral. 



5. A. Mapping the Exact Spectral Integral to the Canonical Form 

The canonical integral is constructed about a a typical reference point along the beam 



trajectory, which is taken to be on the edge of caustic 2 (see Fig. 10 ; recall from the discussion 
in Sec. [3] that the caustic topology near the beam axis is hyperbolic umbilic, consisting of 
two sheets denoted as caustic 1 and 2, each has an edge in the symmetry plane y = 0, see 
also Figs. |4]-[5]). This reference point is parameterized, conveniently, by 6*^, the angle of the 
tangent to the edge at this point with respect to the z-axis. The coordinates at this point 



^Specifically, in the present case, one applies a change of variables such that the stationary points of the 



phase function r in the exact spectral integral in Eq. ( 24 ) are mapped to those of the phase <i> in Eq. ( 32 1 , 
while satisfying the constraint that the phases at the respective points are equal, i.e., lot — $. 
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are given by {x2{0z),0, Z2{9z)) of Eq. (A. 15) (see also Eqs. (A. 17) and (A. 19)). The Cartesian 
coordinates around this point, denoted as (x, y, z), are chosen such that z is the tangent to 



the edge and {x,y) are the normal and the bi-normal, respectively (see Fig. 10), and are 
given, therefore, by 



(x, zY = R0^{ 



2 1 ) 



y = y, 



(33) 



where Re is the rotation matrix in Eq. ([t]) with ip = 0^- 

The spectral integral in the (x,y, z) system is obtained from the one in Eq. (24). We start 



by expressing the latter in the (x, y) system of Eq. (|6| by rotating the spectral variables 
such that (f, = R45o(^, with R given in Eq. (|7|). Substituting into Eqs. (|24|)-(|25|) 
we obtain 



2/3 



(27r)^ 



where 



r(f , ?7) = + i&f/?, + (f + i&)f]/c + Ix + fiy + (z, 
A = 0/2c)'/\ 

with ^ = (1 — — f]^)^/^. Recall that here we assume a = 0. 

The spectral variables (^, 57,C) corresponding to {x,y,z) are related to {i,fiX) via 

{icy = ReScY, V = V, 



(34) 



(35) 
(36) 



(37) 



with C = (1 — — ?7^)^/^. Substituting Eqs. ( [33| and (37) into (35), expanding ( to a Taylor 
series about (^, ?)) = (0,0), and keeping terms up to third order, we obtain the following 
expression for r 



T = Ta + (3D,iyc + (3Dilf/c + ix/c + f]y/c + z/c - f z/2c - f (z - 5) /2c, (3^ 



where 



Ta = 2/3sin^^(l -4sin2 ej3)/c, 
A = cos^^(5cos^e^ -3)/6, 



S = 2(3sm^e„ (39) 
= cos9,{l + cos^9,)/2. (40) 



Ta in Eq. (39) is the spectral delay at ;z = for (^,17) = 0. Referring to Fig. 



10 



5 is identified 



as the distance along the tangent from the reference point (x, y,z) = to the cusped edge of 



caustic 1. From Eq. (37) in Eq. (34) is expressed as d^ = (cos 6*2 — sin 6*^ ^/C) d^ ~ cosO^d^ 



where the last term is the approximation for ^ = 0. For ^ 1 one may use 
and D^^l- 91 



el 
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Concentrating only on the z = plane perpendicular to the edge, and normalizing the 



coordinates, the integral in Eq. (34) with the phase r of Eq. (38) can be written in the form 



U = Ae''^^'' U{x,y; S) cos 9, 



where Ta is given in Eq. (39) and U is the canonical integral in Eq. (32) with 



Co = 3i/62^/V[Ai 



The normalized spatial and spectral coordinates in Eq. (|42j) are 



and the normalized parameters are 

5 = {k(3y/^Dl/^D;^6/2(3 = {kpY^^Df^D;^ sin^ 6, 



{k(3/?,f/^& 



(41) 

(42) 
(43) 

(44) 
(45) 

(46) 
(47) 



where in the second equality in 5 we substituted 5 from Eq. (39), while the last term applies 
for 9-^ < 1. 



The spectral phase $ in Eq. (42), which is the simplest polynomial that describes the AiB 



in the vicinity of the caustic, is of the hyperbolic umbilic type 24, pp. 66], |28| pp. 41]. It 
follows that the field-scale near the caustic is governed by the factor (fc/3)^/^ that scales the 
local coordinates {x/(3^y/(3) in Eq. (44). Note that the fc^^^ scaling is expected in view of 



the general theory of diffraction catastrophes 29, Sec. 36.6] 



The caustic topology is controlled by the parameter 5 of Eq. (46) which represents the 
normalized distance between caustics 1 and 2 such that it scales with frequency like {k(3Y^^ 
and also grows with the observation range z. For 5^1, caustics 1 and 2 are far apart and 
may be described by separated cusp and fold catastrophes, respectively. In this case, the hy- 



perbolic umbilic integral may be reduced to the Pearcey and Airy integrals 29, Sec. 36.2(ii)] 
that represent, respectively, the field near these catastrophes. 



Equation (41) describes the AiB field in a plane normal to the beam axis. The solution in 



the entire domain along the axis is obtained by tracking the local (£, z) coordinate system 



of Eq. (33) along the beam axis. The parameters in Eqs. (39)-(40) are therefore expressed 
[z) which describes points on the beam trajectory (see Eq. (A. 15)). Since this 



in terms of 9z 

local spectrum representation is constructed entirely using the ray system topology near the 
caustic, it follows that the ray approach can be used to construct AiB-type solutions that 
follow any prescribed convex trajectory by back-projecting the rays to construct the relevant 
aperture distribution (e.g., the 2D example in [15j). 
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There are several alternative forms for the hyperbolic umbilic spectral phase all sharing 
the "germ", i.e., the highest order terms of $ which depend only on {C,,'f]) and define the 



catastrophe. The phase $ in Eq. (42) contains the germ ^ + ^rj , but there is another form 



used in 29, Sec. 36.3] which is discussed in Appendix Bl The rest of the phase terms, called 



perturbation terms (or unfolding terms), can have different forms 35, chap. 4] without 
changing the catastrophe type, as long as they are small. For example, the phase used 
in 28, pp. 41] is the same as in Eq. (42) with f/^(5 replaced by ^^5, but since these are 



perturbation terms, the two forms are equivalent. In fact, the latter can be obtained by 
setting the reference point of the local coordinates on the edge of caustic 1 and taking the 
z axis along the direction of the ray that touches the cusped edge there (see point z = 6 in 



Fig. 10). We prefer the form in Eq. (42) since in the alternative one, z is not tangent to the 



edge of caustic 1 (see dashed line in Fig. 10) and the tracking along this edge becomes more 



complicated. Finally, in Sec. 5.B below we extend $ to the finite energy AiB case. 



5.B. Finite Energy AiB 



Referring to Eq. (35), we now consider the case where a 7^ 0. Repeating the procedure 



thereafter we obtain a polynomial approximation for r which is similar to Eq. (38) 



t = t: + PD^ e/c + PD^ /c + ^xyc + vy/c + z/c -ei^- e)/2c - fj^z - r)/2c. 



(48) 



except that the coordinate x of Eq. (33) and the coefficients r^, D12 and S of Eqs. (39)-(40) 
are now replaced by 



x" = X + sin(26'^) — cos6'2], 
r " = + f3[ia sin^ 9^ — o? sin 6*2 — ia^/3]/c, 
5" = 5 + /3[za; + ia cos^ + 0? sin 6*2], 
= Di_2 — ia sin Qz cos d^j c. 



and the term ^z in Eq. (38) now becomes t^[z — e) with 



e = /3(za + a sin 6^2/2 — sin 6^/2). 



For z = 0, $ of Eq. (42) becomes 



$(e, r?; X, y- S; e) = f + ^t]^ + ^x + rjy + f e + ri'^S 



2- I -21 



(49) 
(50) 
(51) 
(52) 



(53) 



(54) 



where 



(55) 
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and the other normahzed variables are given by Eqs. (45)-(46) with x, S, Ta, D^ 2 replaced by 



the corresponding parameters x", 5", r°, D'^2 of Eqs. @-rt52|. Note that Eq. d54| reduces 
back to Eq. (|42| for a = 0. 



Expression (41) for the field becomes 



U = Ae'""^" U{x,y] 6; e) cos 9, 



(56) 



where U is the canonical integral in Eq. (32) with $(^, rj; x, y; 6; e) given in Eq. (54). 

The expression in Eq. (54) has an additional perturbation term ^^e (assuming a is small ), 



but as noted at the end of Sec. 5. A the diffraction pattern of the field will remain hyperbolic 
umbilic, as can be verified both theoretically and experimentally fTl|5|. The main effect of a 



is an overall amplitude decay which is mainly due to the term e*'^^° in Eq. (56). 



5. C. Numerical Evaluation of the Canonical Integral 



The canonical integral in Eq. (32) with $ in Eq. (42) or (54) can be evaluated by a straight 



forward numerical integration. To accelerate the convergence, the integration contours for 
1^1, \ri\ ^ 1 can be deformed into the complex plane so as to reach infinity along the asymp- 
totic valleys of exp(i$). Fastest convergence is obtained if the ^-integration starts asymp- 
totically along the line arg,^ = 57r/6 and ends along the line arg^ = 7i/6. Likewise, the r/ 
integration may start asymptotically along the line argr/ = — 37r/4 and ends along the line 
argf] = tt/A. Another approach is to use the ID integral representation for the hyperbolic 



umbihc as explained in 29, Sec. 36.15(iii)]. 



Here we shall use the standard series expansion for the canonical double integral in 29 



Eq. 36.8.3] and extend it to the case of finite energy AiB. The small parameters in the series 
expansion are S ~ 0{{k$Y^^ sin^ 9 z) and e ~ 0{{k$Y/^a). For typical values of a,u, $ such 
that 5, e <^ 1 we may expand into a Taylor series the term exp[i{C,'^e + fj'^S)] appearing in the 



expression for U in Eq. (32) with $ given in Eq. (54). We obtain 



U{x,y;6; e) 



p,m 



\ 2m i^{^,ri;x ,y) 



(57) 



where 



^{L V, X, y) = $(^, ri;x,y;5;e) - ^ e-r] S 
= P + ^r/^ + ^x + f]y. 



(58) 



Using {iS,) ^ dx.iirf) <(-)■ dy in Eq. (57) and taking the derivatives outside of the integral, we 
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rewrite Eq. (57) as 



(59) 



I{x,y) 



d^df] e 



with {u,v) given by 



(27r)2 

{u,vy = {3/2Y/'Rl,o{3-'/'x, y)- 



(60) 



(61) 



where R is given in Eq. 0. Equations (|59|-(|6T| can readily be calculated using standard 



numerical routines. Note that in the limit e — )■ 0, the series representation in Eq. (59) reduces 



to the standard series expansion in 29, Eq. 36.8.3] 



5.D. Numerical Example 
and 
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depict the field intensity at two different frequencies such that k(3 



Figures 

10^ and 10^, respectively. The results are shown along the normalized x-axis in the y = 
and y = 1 planes. The observation domain is centered at the point {x,z) = (0.006, 0.16)/3 
on the edge of caustic 2. The figures compare the field calculated via the canonical integral 
in Eq. (41) (green line), with the PE solution of Eqs. ([2])-([3]) (blue points), and the GO 



field of Eq. (27) (red dashed line). The canonical integral has been calculated via the series 
expansion in Eqs. (59)-([60|), keeping only the terms p = 0, ...3 and m = 0,...3. A decay 
parameter a 



10 ^ has been used. 



In Fig. 11 the canonical integral solution and the PE solution are in an excellent agreement 
on both the lit and the shadow sides of the caustic. As expected, the GO solution breaks 
on the caustic, but otherwise it agrees well with the other solutions. Further away on the lit 
side of the caustic, the GO solution starts to deviate, but in this region it is, in fact, more 
accurate than the other solutions since it does not rely on the local spectrum approximation 
or on the PE approximation. 



In Fig. 12, on the other hand, on observes a significant translational error of the PE 



solution (notice, on the other hand, the perfect agreement between the canonical integral 
and the GO solutions, except on the caustic). The error of the PE solution is due to the 



deviation of the PE-based beam trajectory in Eq. (10) from the true ray-based trajectory, 



and can be parameterized by Xp which is the normalized location of the PE caustic in the 



x-system (recall from Fig. 10 that the x-system is centered on the edge of caustic 2, see also 
Fig. |4]). Xp can be found using Eq. (10) in conjunction with Eqs. (33), (A.14)-(A.15) and 



(44). Here we consider only the small 9z approximations of Xp, obtaining 



Xp ~ 



(62) 
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In Fig. 11 the frequency is relatively low such that Xp = —2.7 x 10 ^ (i.e., Xp = —6 x 10 ^A), 



hence the error of the PE solution is hardly noticeable on the figure scale. In Fig. [I2| on the 
other hand, the frequency is higher giving Xp = —0.59 (i.e., Xp = — 6A). Since Xp is 0(1), the 
translational error of the PE solution is readily discerned in the figure. 

In addition to the translational error of the PE solution, the exact canonical integral has 



some intricate features which are readily discerned in the main-lobe of Fig. 12 (notice the 



peak 1% error of the PE-solution there). Recalling the discussion after Eq. (46), the field 
structure of the canonical integral is parameterized by 6 which is the normalized distance 
between caustics 1 and 2. As noted there, this parameter grows with both frequency and 
observation range. For 6 <^ 1 (near the aperture or small k(3), the distance between these 
caustics is small and the canonical integral reduced to the Ai function field structure in 
Eqs. ([2])-(|3|, but for larger 5, the field structure of the canonical integral describes the 
separation between the two caustics and is more complicated than the PE solution. This is 
seen in Figs. 11 and 12 where 6 = 7.6 x 10~^ and 3.5 x 10~^, respectively. In these figures, the 
cusped edge of caustic 1 is located at Xi = —1.7 x 10~^ and Xi = —3.7 x 10~^, respectively 
(i.e., xi = -4 X 10-^A and = -4 x 10^^^). 

Finally, to help the reader we map the parameters used in these numerical examples to the 
conventional ones used for example in [2|[5] (see also discussion after Eq. ([T]) ) . Noting that the 
"Fresnel distance" Zn 



kxl used there is related to our P parameter via kz^ 

is z 



it follows that the range at the caustics in Figs. [TT] and [12 
respectively. 



2^l\k~Pfl\ 
2.7^0 and z = 12.5^0, 



6. Discussion on the Complete Caustic Topology 

We consider the global caustic topology which includes far away parts of the caustic which 



were not included in the analysis of Sec. 3.B The field in these regions is formed by contribu- 
tions of rays that emanate sideways from far away points in the aperture. This global picture 



is presented here in connection with the discussion in Sec. 6.C regarding the limitations of 
the AiB. Since the caustic structure is rather complex, we discuss separately the topologies 



of caustic 2 and 1 in Sees. 6. A and 6.B, respectively. 



6. A. Global Structure of Caustic 2 

The global structure of caustic 2 is depicted in Fig. 13 , as a zoom-out view of Fig. |5| In 
addition to the sheet shown in Fig. [5} tagged here "caustics 2A" , there is another caustic 
surface which is tagged "caustic 2B." These two surfaces are connected and form a pyramid- 
like structure which terminates at a sharp edge at the top. Caustic 2B is formed by rays 



exiting the aperture at far away regions, as illustrated in Fig. 16 which shows a cross section 



of the caustic surface in the y = (the cross sectional line is also depicted in Fig. 13 by a 
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red dashed-line) . 



6.B. Global Structure of Caustic 1 



Figs. 14(a) and 14(b) depict the complete structure of caustic 1 as a zoom-out view of Fig. [5| 
with Fig. |14(a) showing the caustic parts formed by ray species s = 1 (see Fig. [2] for the 
species definition), and Fig. |14(b)] showing the rest of the caustic. The surfaces in Figs. 14(a) 
and [14(b) are connected at the margins that are dehneated in the figures. The connections 



of these surfaces are clarified in Fig. 15 which depicts cross-sectional views of the complete 
caustic at two constant-2; planes. A cross sectional view in the symmetry plane y = is 



depicted in Fig. 16 



As can be discerned from Fig. 14(a) , the part of caustic 1 that is shown in Fig. |5| tagged 
here as "caustics lA", is connected to a hat-like structure, tagged as "caustic IB." The latter 
is formed by rays of species 1 exiting the aperture from far away points, as in the case of 



caustic 2B of Fig. 13 Caustic IB is a cuspoid (as may also be seen in Fig. 15 ) with its edge 



in the y = plane (red dashed line in Fig. 14(a)[ ) 



6. C. Discussion: Limitations of the AiB Solution 



As discussed in Sec. 3.B the AiB has been identified as a hyperbolic umbilic diffraction 



catastrophe. However, as can be seen from Figs. 13 14(a) , near the termination point of 



caustic 2, the hyperbolic umbilic caustic is no longer isolated, and the caustic structure 



evolves into the more complex parabolic umbilic catastrophe 24, Fig. 7.3]. From this global 
analysis it then follows that the propagation range of the AiB is limited to the termination 



point of caustics 2A (see also Fig. 16) where the AiB loses its beam shape. The maximal 
propagation range is of order (5 and hence it can be controlled by choosing /3 for a desired 
application. The range is also controlled by a which controls the decay rate along the beam 



axis (see Eqs. (56) and (50)) 



The exact ray-based solutions of Sees. |4] and |5] may be used up to the maximal propagation 
range discussed above. In the near zone, one may also use the PE solution in Eqs. 
but the accuracy of this solution reduces as the distance between the exact and the PE 
approximated caustics grows (see Fig. [6]), giving rise to the translational error of the PE 



field in Fig. 12 



6.D. Comments on the Range Limitations of the 2D AiB 



In 11 and 16 it has been shown that the 2D-AiB is supported by a smooth caustic which 



is a fold catastrophe 28 . Those references, however, have concentrated on the field near the 



beam axis. A global analysis reveals that this 2D caustic is actually a cusp whose form is very 



similar to the cross section of caustic 2 in Fig. 16 (caustic 1 does not exist in the 2D case 
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This global structure was not discussed in 11 and 16 since the rays forming the missing 



parts have weak contributions and late arrival times beyond the pertinent time window 16 



The analysis in Sec. |6.C| regarding the limitations of the AiB can be similarly applied to the 
2D case. 



7. Conclusions 

We presented a wave-theoretic analysis of the 3D- AiB and derived uniform solutions for the 
transition layer near the propagation axis. The formulation is based on exact spectral and 
ray representations and therefore does not suffer from the limitations of the conventional 
PE approximation. As in the 2D-AiB case, the 3D-AiB is described by sideways radiating 
rays that focus onto a caustic that delineates the beam trajectory. However, the 3D ray 
topology is considerably more complicated than the 2D topology, and it has many intricate 
features. Specifically, near the curved beam axis, the ray topology of the 3D-AiB consists of 
a cusped double-layered caustic identified as the hyperbolic umbilic catastrophe (Sec. |3|. At 
a certain range away from the aperture, the hyperbolic umbilic caustic is no longer isolated 
and evolves into the more complex parabolic umbilic caustic, where the AiB loses its beam 
shape. 

To clarify these issues we have investigated the ray topology of the 3D-AiB using the 
exact spectral representation for the AiB. We concentrated first on the region near the 
curved beam axis (Sec. |3]), but we have also explored the global caustic topology (Sec. [6]). 
Specifically, we have thoroughly explored all the ray paths that form the double-layered 
hyperbolic umbilic caustic noted above (Figs. [7 10). These contributions are essential in 
order to understand the time-dependent Airy pulsed beam (AiPB) that will be presented in 



a subsequent publication 21 



The field calculations in Sees. |4]-[5] were based on an asymptotic evaluation of the exact 
spectral integral of the AiB. In order to calibrate the final results, we derived first the lowest 



order geometrical optics (GO) solution (Sec. 4.B) which, however, fails near the beam axis. 
The field near the beam axis has been considered separately in Sec. |5] within the framework 



of catastrophe theory 24,27,28 . As discussed in the introduction (see also Sec. pi), a globally 



uniform field solution can be obtained via the procedure in 30, Sec. 5]. This approach, while 
very elegant and systematic, leads to implicit representations where the field is described 
in a generalized coordinate system that has to be calculated numerically. We therefore used 



in Sec. 5. A an alterative local approach where the spectral integral of the AiB has been 
approximated near the caustic using a local coordinate system, and the final result has been 
expressed as a canonical hyperbolic-umbilic diffraction integral [29^, Sec. 36.2]. This more 



direct approach does not have the global uniformity of the general implicit procedure in 30 



yet it has led to a solution that is expressed explicitly in the local coordinate system of the 
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caustic, where all the field parameters are expressed in terms of the geometrical parameters 
of the caustic and of the aperture field distribution. This solution is valid in the transition 
layer near the caustic and reduces correctly to the GO field away from the caustic, but far 



away, one has to switch to the GO field of Sec. |4.B| noted above. Using this approach we 
also modified the hyperbolic-umbilic integral by introducing a complex loss parameter that 
describes the finite energy AiB. 

The final field solution in Eqs. (56) with (32) and (54)) provides an explicit solution in 
the transition layer near the beam axis. The width of the transition (or boundary) layer 
is of order k~'^^^$^^^ (see x, y coordinates in Eq. (44)) where k is the wave-number and 
/3 is a frequency-independent parameter describing the beam geometry (see Eq. (10)). The 
caustic topology is controlled by the dimensionless parameter 6 of Eq. (46) appearing in 
$ of Eq. (54). This parameter represents the normalized distance between the two layers 
of the hyperbolic umbilic caustic, scaled with the normalized frequency like (/c/3)^/^. As 5 
grows with the observation range z, the structure of the boundary layer field described by 
the canonical integral becomes more complex while for 5 ^ 1, the two layers are separable 
and the field can be described by separate cusp and fold catastrophes, as discussed after 
Eq. (|46|). 

The canonical hyperbolic umbilic integral can be evaluated efficiently via contour deforma- 



tion in the complex plane. Here, however, we calculated it using the series expansion of 29 



Sec. 36.8] which has been modified to accommodate the finite energy AiB (Sec. |5.C ). We 
used it to calculate the field near the beam axis under several parameter regimes (Sec. 5.D). 



The accuracy of the canonical integral has been demonstrated by noting in Fig. [TTJthat away 
from the beam axis it reduces to the GO field (the GO field near the hyperbolic umbilic caus- 
tic is a sum of four rays). Far away from the transition layer, however, the field is described 
more accurately by the GO solution of Sec. 4.B since it is structured upon the exact ray 
skeleton and not on the local ray skeleton near the beam axis used here for the construction 
of the transition function (note that the solutions obtained via the globally uniform approach 
in 30, Sec. 5] reduce correctly to the GO field even very far from the transition layer, yet. 



as noted earlier, these solutions are implicit functions of the space coordinates). Comparing 
the canonical integral to the conventional PE solution one observes in Fig. [12] a significant 
translational error of the latter which is due to the error of the PE approximated beam 
trajectory. This error is parameterized by the dimensionless parameter Xp of Eq. ( 62 ) which 
describes the distance between the exact and the PE caustics. Xp grows with z (see Fig s. [4| 
andp| and scales with the normalized frequency like (fc/3)^/^. As discerned from Figs. 



Fig 


s.p 
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the error is small for Xp ^ 1, but it becomes significant for Xp ~ 1. 

In summary, it has been demonstrated that the local canonical integral in conjunction 
with GO provides an accurate and convenient solution to the AiB, as an alternative to the 
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conventional PE solution. This solution is structured upon the exact ray skeleton and there- 
fore has a wider validity range than the PE solution. Furthermore, this solution provides 
a convenient design tool for synthesizing other types of beams propagating along any pre- 
scribed convex trajectory in uniform or non uniform medium. An application of this idea in 



2D configurations has been presented in 15 



The analysis in this work has been performed within the framework of the non- dispersive 
AiB where the aperture field is scaled with frequency such that the ray skeleton is frequency- 
independent. We use the non- dispersive formulation because it simplifies the ray analysis 
by making it more transparent geometrically. Moreover, this formulation also extends the 
AiB solution to the ultra wide band (UWB) frequency regime by ensuring that the pulsed 
field propagates along the curved beam trajectory without dispersion. This non-dispersive 



formulation will be used in a subsequent publication 21 to introduce a new class of 3D Airy 
pulsed beams (AiPB). 

A. Appendix A: Derivation of the GO Field from the Spectral Representation 

In Sec. [3] we introduced the ray representation of the AiB based on the approximate initial 



field distribution in Eq. (15). In this appendix we derive the ray representation based on the 



exact spectral representation of the initial field distribution in Eq. (24), and show that the 
results are the same as those of Sec. El 



As discussed in Eq. (15), a is relatively small and hence we consider first the ray and 



caustic equations for a = (Sees. A.1-A.3) and then we consider the GO field for the case 



a^O (Sec. A.4). 



A.l. The Ray Equation 



We start by deriving the ray equation in Eq. (11) from the saddle point condition of the 



exact spectral representation in Eq. (24). 



Substituting the spectral delay function r of Eq. (25) with a = into the saddle points 



condition in Eq. (26) we obtain 



% = %o + x/c-(^/c)(e/C) = 0, 
dr,T = dr^To + y/c - {z/c){r]/C) = 0, 



where ( is defined in Eq. (13). For z = 0, Eqs. (A.1)-(A.2) reduce to 



dr/To 



-x'/c = 13^"^ /c, 
-y'/c = l3r]'^/c, 



(A.l) 
(A.2) 



(A.3) 
(A.4) 



where {x',y') are the coordinates in the z = plan. The right hand side in Eqs. (A.3) 



(A.4) is obtained by using Eq. (22) for Tq and it is identical to Eq. (18) for the ray exit 
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angles. Substituting the left hand side of Eqs. (A.3)-(A.4) back into Eq. (A.1)-(A.2) and 
rearranging, we obtain 



X = x' + ^{z/C), 

y = y' + r]{z/C), 



(A.5) 
(A.6) 



which are, in fact the expressions for the ray trajectories in Eq. (11) with a = z/(. It follows 



that the ray skeleton derived in Sec. |3j using the approximation in Eq. ( 14 ) is the same as 
the one derived here from the exact spectral representation. 

A. 2. Ray and Caustic Topology - Exact Analysis 

There are two alternative routs to derive the expression for the caustic of Sec. |3| The first 



is to calculate JT" of Eq. (19) using the ray trajectories in Eq. (11), and then to solve for its 



zeros. The second relies on the spectral formulation of Sec. 4. A Here we start with the latter 
since it also complies with the formalism of catastrophe theory which is discussed in Sec. [5j 
To calculate the caustic surface we note that the elements of the Hessian matrix H in 



Eq. (28) are given by 



(c//3)a|r = 2e-(^//3)(l-r^^)/C^ 
(c//3)9jr = 2r/-(^//3)(l-a/C', 



and hence, the determinant H is given by 



H =[Ciz' + C,(3z + C,(3^]/c\ 



(A.7) 
(A.8) 
(A.9) 



(A.IO) 
(A.ll) 
(A.12) 



At the caustic, two or more rays (saddle points) coalesce, implying that H = there. This 



equation should be solved in conjunction with Eqs. (A.1)-(A.4) since the caustic is composed 



of rays. From Eq. ( |A.10[ ) a ray with parameters (^, 77) is tangent to the caustic at two points 

(A.13) 



The {x,y) coordinates of these points are found by substituting Eq. (A.13) back into 



Eqs. (A.1)-(A.2). 



As shown in Sec. [3] and Fig. [5} the exact caustic topology near the beam axis is composed 
of two sheets, denotes as caustic 1 and caustic 2. We have also shown in Figs. [7]-[8]that a ray of 
species 1 touches caustic 2 first and then touches caustic 1, never touching the caustic again. 



Thus the "±" sign in Eq. (A.13) corresponds to a point on caustic 1 and 2, respectively. 
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The edges of the caustics in the symmetry plane y = (see Figs, [s] and 111) are found by 



setting ^ = ?7 in Eqs. (A. 13) and (A.1)-(A.2). The {x,z) coordinates can be expressed in the 
following parametric form 



(a;i,2i) = /3(sin^6'^, sin 26*2), 

(xa, Z2) = /3(sin^ 9z cos 26*2, sin 26*2 cos^ 9z 



(A. 14) 
(A.15) 



where sin^^^ = (,^^ + rj'^Y^^, 9^ being the ray-angle with respect to the z axis (see Eq. (12)). 
Alternatively, 2:12 in Eqs. (A.14)-(A.15) can be expressed as function of x 

Vx/Pf'')/^ 



zj(3 = 2{i/(3fl\l-~x/(3fl\ 



(A.16) 

z,/~P = 2q^'\l - qfl\ g = (1 - (1 - 8£/^)^/2)/4, (A.17) 

where all square roots taken with a positive real part. At moderate ranges such that x/ (3 
these expressions may be approximated by 

^ 2{i/~pfl^{l - i/A~P), (A.18) 
zS ^2{i/~pf'^{l-2>x/2~P), (A.19) 

where we used (1 



X] 



! 1 — nx. The lowest order approximation in Eqs. (A.18)-(A.19) is 
the PE solution Zp/13 ^ 2(i//3)i/2 of Eq. (|lO|), hence one finds for a given x that Z2 < z^ < Zp 
(see also Fig. [6| . 



Finally we note that the same results could have been obtained via the ray-system analysis. 
The Jacobian in Eq. (20) can be expressed explicitly as 

J{a) = CC'a^ + PC^a + (5^C,C, (A.20) 

where Ci_2,3 are given in Eqs. ( |A.ll )-(A.12). Comparing Eqs. (A. 10) and (A.20), it follows 
that the zeros (Jci,2 of J {a) correspond to the zeros Zcr,2 of H in Eq. (A. 13). (Tci,2 = Zcx-ijC, 
define the ray lengths to the caustic. 

A. 3. The Caustic within the Parabolic Equation Approximation 

In the paraxial regime one can approximate in Eq. (13) C ~ 1 ~ ^^/2 — ^7^/2. Substituting 
C ^ 1 in the ray equations (|A.ll)-(|A.2l) with (|A.3l)-(|A.4|), yields 



e = z/2f3 ± V{z/2f3)^ - x/f3, 
r, = z/2f5 ± ^{zl2f5Y-ylp. 



The solutions of Eq. (A. 10) yields two surfaces S'i_2, given by 

S,:{z/2fif-xll3 = Q & y/l3>{z/2l3)\ 
S2:{z/2(3f-y/(3 = Q & x / (3 > {z /2(3)\ 



(A.21) 
(A.22) 

(A.23) 
(A.24) 



which are shown in Fig. [3j These surfaces intersect at the edge which is the beam axis in 
Eq. (|T0|. 
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A. 4- Details of the Derivation of the GO Field in Eqs. (27)-(31) 



The phase term r,. in Eq. (30) is obtained by substituting the ray directions in Eqs. (A. 3) 



(A. 4) (see also Eq. (18)) and the ray equation in Eq. (12) into t(^,?7) of Eq. (25). 



Next, in order to derive the amphtude term in Eq. (31), we note from Eq. (A. 10) and 



Eq. (A. 20) that the determinant of the Hessian and the Jacobian can be cast in the form 



Hia) 



(A.25) 



where a = z/C^ and aci,2 correspond to Zci,2 of Eq. (A. 13); crci,2 are the ray lengths to caustic 1 
and 2, respectively. 



From Eq. (A.25), the eigenvalues of the Hessian matrix defined in Eq. (29) are given by 



hi,2 = (o- - o-ci,2)/c. 



(A.26) 



Thus /ii < if the ray has not touched caustic 1 yet, while /ii > if it has touched it. 
Similarly, /12 ^ if the ray has not yet touched or has touched caustic 2, respectively. We 



also note, by analyzing Eq. (A. 13) for 2ci,2 using Ci_2,3 from Eqs. (A.ll )-(A.12), that the ray 



species s = 1 has (Tci,2 > 0; species s = 2, 3 have < and > 0; and species s = 4 has 



c"ci,2 < 0. Combining these results we find that of Eq. (29) can be written as 



/ir = sgn{/iir} + sgn{/i2r} = 2(-l + + 



(A.27) 



where M^g is specified after Eq. (17) and is the Maslov index defined after Eq. (31) that 



counts the number of times that the rth ray has touched a caustic. 



Next, using the left hand side in Eq. (A.25) in conjunction with Eqs. (A.3)-(A.4) we find 



that in Eq. ^ can be expressed as 



\H\-^" = c2-'r^'\-xr"\-y'V\J{^)IJ{<y)\"' 



(A.28) 



The final result for the GO amplitude in the right hand side of Eq. (31) is obtained by 



substituting Eqs. (A.27)-(A.28) in the left hand side of that equation. 



A. 5. Finite Energy AiB (a ^0) 



For a ^ the saddle points {^,ri) are complex and are given by Eqs. (A.1)-(A.4) with ^ 



and 7] in Eqs. (A.3)-(A.4) replaced hy ^ + ia and r] + ia, respectively. Since a is small, the 



solution for the complex saddle points is very close to the real saddle points {^r, Vr) obtained 



for a = 0. Therefore, we expand ( in Eq. (A.1)-(A.2) to first order about {^r,Vr 



(A.29) 
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where Cr = C; C| = d^Cy Q = "^^C; all taken at rj) = (^j,, r^,.), and are given by 

0.=-^r/Cr, C = -WC.. (A.30) 



Substituting into Eqs. (A.1)-(A.2) and solving for the saddle points {^,1]) yields 



(A.3i: 



Substituting Eq. (AJU), and {AJ29)-{AM) into r of Eq. (25), using Eqs. (11) and (A.3)- 



(A.4) and neglecting terms 0{a ) and lower, we obtain 



(A.32) 



Tr = CTor + car IS the expressiou in Eq. (30). The term + in Eq. (A.32) appears inside 



the exponent in Eq. (17) where it is interpreted as a factor in the initial amplitude A^s rather 



then a complex time delay. In the expression for the Hessian we neglect a so Eq. (A. 28) is 
without change. 

B. Appendix B: Additional Canonical Form for the Hyperbolic Umbilic Catas- 
trophe 

In order to prevent possible confusion we note that the Hyperbolic Umbilic catastrophe 



has another canonical form in addition to the one in Eq. (|42|), which has the following 

(B.l) 



form [29| Sec. 36.3, Eq. (3)], (Mj pp. 37] 

v; X, y; 6) = f + f + ^7)6 + + T)y. 
This form is obtained from Eq. (|42|) by using the following transformations 



(x,y)" = 2-/«R4V([x + 35V4],3^/^y)^ 
5 = -36/2^/\ 



(B.2) 
(B.3) 
(B.4) 



where (x, y) and 6 are defined in Eqs. (44)-(46), (^, rj) are defined in Eq. (45), and R is given 
m Eq. 0. Note that some terms which do not depend on (^, 77) have been omitted since 
they do not affect the catastrophe type. 
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Fig. 1. (Color online) Field intensity |f/p at (a) the aperture plane z = 0, 
and (b) the z = 0.8/3 plane. The field is given by the PE solution in Eqs. ([2])- 
(|3|. Beam parameters are as in (Q-Q with kf3 = 100, a = 10"^ The beam 
propagates in the plane of symmetry y = in the direction = 45°. The two 
coordinate systems {x, z) and (x, y) are related via Eq. Q. The zeros = 
in Fig. (a) along the x and y axes are given, respectively, by x = an/3{(3k)~'^^^, 
y = a„/3(/3A;)~^/^, where with n = 1,2,..., are the zeros of Ai given in 
Table 9.9.1]. 
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(c) Ray species s = 3 



(d) Ray species s = 4 



Fig. 2. (Color online) Projections of the ray directions s defined in Eqs. (12) 
and (18) onto the z = plane (white arrows). Also shown for reference the 
intensity |f/p from Fig. 1(a) Figures (a)-(d) show the four ray species defined 
in Eq. ^ such that (a) 6^, By G [0,7r/2]; (b) 6^ G [7r/2,7r] and By G [0,7r/2] 
(c) 9, G [0,7r/2] and Oy G [n/2,n] (d) 9^,9y G [n/2,7i]. 
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Fig. 3. (Color online) The caustic of the AiB within the PE approximation 
(gray surface). The caustic is analyzed in Appendix A. 3, It consists of two 
curved sheets S'i_2, defined respectively, in Eqs. (A.23)-(A.24) (dark and lit 
sheets in Fig. (a), respectively). The sheets are joined at the edge (dashed red 
line) which delineates the beam propagation trajectory. Cross-sectional cuts 
are shown in Fig. |4} 
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Fig. 4. (Color online) Cross sections of the caustics of the AiB at two different 
ranges: (a) z/P = 0.2; and (b) zjP = 0.25. Black line: the caustic of the 
PE approximation in Fig. |3j Gray and green lines: the exact caustic, which 
consists of two surfaces denoted as "caustic 1" and "caustic 2" respectively 



see Sec. 3.B.2 and Fig. k^. The x axis in this figure is centered about Xi{z) 



the edge of caustic 1 in Eq. ( A.16 ). The insets zoom on the area near the origin 
where caustic 1 has a "cusped edge" while caustic 2 has a smooth corner. One 
observes that the distance between caustics 1 and 2 and the PE approximated 
caustic increases with z (see also Fig. [6]). 
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Fig. 5. (Color online) The exact caustic of the AiB, discussed in Sec. 3.B.2 



and in Appendix |A.2 It consists of two surfaces denoted as "caustic 1" and 
"caustic 2" (gray and green surfaces, respectively). The distance between these 
surfaces is hardly noticeable on the scale of this Figure, but it is shown in the 
cross-sectional cut of the caustics in Fig. |4j Each caustic consists of two surfaces 
joined at an edge (dashed red lines). The ray formation of caustics 1 and 2 are 
described in Figs. [7] and |8| respectively. 
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Fig. 6. (Color online) The difference between the edges of the PE approximated 
and of the exact caustics. The figure plots Zp — Zi as a function of x, where Zp 
is the coordinate of the paraxially approximated beam trajectory in Eq. (10), 
and Zi, i = 1,2, are the coordinates of the exact edges of caustic i, given in 



Eqs. (A.16)-( A.17). From Eq. (10 ), the horizontal range x/ j3 = 0.1 corresponds 
to the observation range Zp/$ = 2\/0A and there {zp — Zi)/zp ^ 5% — 6%. 
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Fig. 7. (Color online) Ray formation of caustic 1 (gray surface) of Fig. [5] 
Each exit point (x', y') in the aperture (red point) emits four rays (black lines) 
belonging to species s = 1,...4 of Fig. |2} Ray 4 does not touch the caustic. 
Rays 2,3 touch only caustic 1 (red points). Ray 1 penetrates caustic 1 (through 
the white circle), touches caustic 2 (green point), penetrates caustic 1 again 
(through the second white circle) and finally touches caustic 1 (red point), 
never touching the caustic again. The figure shows traces of the exit points 
in the aperture (thick lines) and the corresponding traces of the points of 
tangency on the caustic (thin lines with the arrows indicating the direction of 
the traces). To emphasize the symmetry, traces with y' ^ appear as full or 
dashed lines. Lines of constant y' and of constant x' are plotted in magenta 
and blue, respectively. The traces corresponding to the edge of caustic 1 are 
shown as red dashed-lines. The parameters used here are: x[ = —0.004/3, x'^ = 
-0.009/3, y[ = ^0.0015/3, y'^ = ^0.0065/3, y', = 0. 
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Fig. 8. (Color online) Ray formation of caustics 2 (green surface) of Fig. [5] 
The notations are the same as in Fig. [7} Note that caustic 2 is formed only by 
ray species 1. 
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Fig. 9. (Color online) Rays reaching an observation point at r = {x,y,z) = 
(0.005, 0, 0.17)/3 (blue point) on the lit side of the caustic, tagged by the ray in- 
dex r = 1, ...4. The parts of the rays before and after r are depicted as sohd and 
dashed lines, respectively, and the parts between caustic 1 and 2 are depicted 
as dotted lines. The ray initiation points {x'^,y'^)/(3 are: ray 1 (—0.14,0)10"^; 
rays 2,3 (—0.95, =F. 81)10"^; ray 4 (—1.8,0)10"^. The ray direction parameters 
{ir.nr) are: ray 1 (2.7,2.7)10-2; ray 2 (2.7,9.4)10-2; ray 3 (9.4, 2.7)10-^; ray 4 
(9.5,9.5)10-2. The paths of rays 2 and 3 are similar to that of ray s = 1 in 
Fig. [7j they penetrate caustic 1 (through the white circles), touch caustic 2 
(green points), penetrate caustic 1 again (upper white circles) and then touch 
caustic 1 (red points) after passing through r (blue point). Rays 1 and 4 prop- 
agate in the y = plane. Ray 4 touches the caustic after passing through r: It 
penetrates caustic 1 (through the white circle), touches the edge of caustic 2 
(green point) and then penetrates back through the cusped-edge of caustic 1 
and at the same point it is also tangent to the cusp (red point; see discussion in 
the last paragraph of Sec. 3.B.2 and also Fig. [Io| ). The path of ray 1 is similar 
to that of ray 4, only it already touched the caustics before passing through r. 
When r moves toward the cusped-edge of caustic 1 from its lit side, rays 1-3 
coalesce. When r is on the shadow side of edges 1 such that it is between the 
edges of caustics 1 and 2, it is reached by 2 rays which penetrate caustic 1, one 
which results from the coalesced rays 1-3 and has already touched caustic 2, 
and ray 4 which has not touched caustic 2 yet. These two rays coalesce as r 
moves toward the edge of caustic 2. 
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Fig. 10. (Color online) Schematic view of the edges of caustics 1 and 2 (dashed 
and solid red hnes, respectively) in the symmetry plane y = 0. The edges 
are given in Eqs. (A.16)-( A.17[ ). As shown in Fig. |4| caustic 1 is actually a 
cuspoid but only its edge is shown here. A ray (dashed black line) penetrates 
through caustic 1 before the point of tangency with caustic 2. After that point 
it penetrates again through caustic 1 (point z = S in the figure). At that 
point it is also tangent to the cusp of caustic 1 (see see discussion in the last 
paragraph of Sec. 3.B.2). Also shown, the local Cartesian coordinate system 
(x, y, z) centered around a point on caustic 2, such that z is the tangent to the 
edge at that point, while x and y = y are the normal and the bi- normal there, 
respectively. This system is parameterized by the angle 6z of the i-axis with 
respect to 2;-axis, and it is used to construct the canonical integral of the field 
in Sec. [5l 
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Fig. 11. (Color online) Intensity |f/p of the field in the planes [a) y = and 
(b) y = I, plotted along x, the normalized of Fig. LLOl which is normal 



to the beam axis, and centered at a point on this axis; here this point is 
{x,z) = /3(0. 006, 0.16). The fields were calculated via the catastrophe theory 
(green line), the PE approximation (blue dots) and GO (red line). The intensity 
is normalized with respect to max |f/op = Ai2(0) of Eq. ^. Parameters: k(3 = 
10^ a = 10"^ 
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Fig. 12. (Color online) Same as Fig. 11 but for kjS = 10' 
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Fig. 13. (Color online) The global structure of caustic 2 (green surface). "Caus- 
tic 2A" is the part that was shown in Figs. [5] and [8} whose edge follows the 
beam axis. This part of the caustic is formed by rays leaving the aperture at 
points close to the caustic. "Caustic 2B" is an additional part of caustic 2 
not shown in Fig. |5] and it is formed by rays leaving the aperture sideways 
from faraway points (see Fig. 16). At some range, caustic 2 terminates at a 
pyramid-like edge. Red dashed-hne: cross section of the caustic at y = which 



is also shown as red dashed-line in Fig. 16 
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(a) The part of caustic 1 formed by species 1 




(b) The part of caustic 1 formed by species 2,3,4 



Fig. 14. (Color online) The global structure of caustic 1. Since the structure 
is complex, we show separately the parts formed by ray species s = 1 and by 
species s = 2,3,4 (Figs, (a) and (b), respectively). In (a), caustic lA is the 
caspoid shown in Figs. [5] and [7] which is formed by rays leaving the aperture 
at points near the beam axis, while caustic IB is formed by rays leaving the 
aperture sideways at faraway points. In (b), caustic IC (blue) and caustic ID 
(magenta) are two overlapping sheets connect to caustic IB at margins 3 and 4 
and to caustic lA at margins 1 and 2 (see also cross-sectional cuts of the caustic 
in Fig. 15). 
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Fig. 15. (Color online) Cross sections of caustic 1 in Fig. 14 at z constant 
planes. 
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Fig. 16. (Color online) Complete cross section of caustic 1 and 2 (gray and 
green lines, respectively) at the plane of symmetry y = 0, also shown as red 
dashed-lines in Figs. 13 and 14 These lines are defined in Eqs. (A.14)-(A.15) 
(or alternatively by Eqs. (A.16)-(A.17)). Black solid lines: rays emanating 
sideways from faraway points, forming caustics IB and 2B. Black dashed line: 
A typical ray near the axis forming caustic lA and 2A (indistinguishable within 
the scale of the figure). Caustics IC and ID of Fig. 14(b) are not formed by 
rays contained in the y = plane and are therefore omitted here. 
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